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Methane storms as a driver of Titan’s dune orientation 
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Titan’s equatorial regions are covered by eastward propagating linear dunes 1 3 i 
This direction is opposite to mean surface winds simulated by Global Climate Mod¬ 
els (GCMs), which are oriented westward at these latitudes, similar to trade winds 
on Earth^l. Different hypotheses have been proposed to address this apparent con¬ 
tradiction, involving Saturn’s gravitational tides®, large scale topography® or wind 
statistics^, but none of them can explain a global eastward dune propagation in 
the equatorial band. Here we analyse the impact of equinoctial tropical methane 
storms developing in the superrotating atmosphere (i.e. the eastward winds at high 
altitude) on Titan’s dune orientation. Using mesoscale simulations of convective 
methane clouds®® with a GCM wind profile featuring super rotation 8 , we show that 
Titan’s storms should produce fast eastward gust fronts above the surface. Such 
gusts dominate the aeolian transport, allowing dunes to extend eastward. This anal¬ 
ysis therefore suggests a coupling between superrotation, tropical methane storms 
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and dune formation on Titan. Furthermore, together with GCM predictions and 
analogies to some terrestrial dune fields, this work provides a general framework 
explaining several major features of Titan’s dunes: linear shape, eastward propa¬ 
gation and poleward divergence, and implies an equatorial origin of Titan’s dune 
sand. 

A major surprise in the exploration of Titan by Cassini was the discovery of large dune fields in 
the equatorial regions, which cover close to 15% of Titan’s surface 13? . These giant dunes are linear 
and parallel to the equator, and are probably composed of hydrocarbon material. The analysis of dune 
morphology around obstacles and dune terminations indicates an eastward dune propagation with 
some regional variations 1 2 (see also Fig. [ 3 J 5 and Supplementary Fig. 5). However, Global Climate 
Models (GCMs) predict that annual mean surface winds are easterlies (westward) at low latitudes 4 , 
as trade winds on Earth®. Therefore, Titan’s dune orientation is opposite to predicted mean winds, 
raising a major enigma. 

Given the non-linear dependence of sediment transport on wind speed, the only way to propagate 
dunes eastward would be the occurrence of episodic fast westerly (eastward) gusts®. Above 5 km, 
Titan’s troposphere is in superrotation with fast eastward winds at any latitude. Pumping momentum 
from this superrotation down to the surface might provide a solution. However, Titan’s tropospheric 
circulation is essentially confined into a 2 km boundary layer 10 . Incidentally, this boundary layer 
circulation may explain the dune spacing of around 2 km 10 u . Dry convection is therefore limited 
to the first 2 km and unable to reach the fast eastward winds 10 . The Koln GCM has been shown to 
generate episodic fast eastward gusts at the equator during Titan’s equinoxes 5 . However, this result 
has not been reproduced by any published GCM, including the Titan IPSL GCM (used here), which 
faithfully reproduces the superrotation and the thermal structure measured by the Huygens probe in 
the troposphere 810 in contrast to the Koln GCM. We therefore conclude that Titan’s atmospheric 
general circulation is not likely to produce fast eastward winds at the equator. 

Similar to the water cycle on Earth, Titan’s weather is characterized by a methane cycle producing 
methane clouds. These clouds are rare over the equatorial region and probably absent during most 
of Titan’s year' 2 . However, during the equinox, the circulation leads to a methane convergence at 
the equator, sufficient to trigger deep convection 1314 . Tropical convective clouds have been detected 
during the equinoctial season 12 15 17 with top altitude between 10 and 30 krrJ 17 ! Huge storm systems 
covering a la rge p art of the equatorial band have been observed, one of them associated with methane 
precipitations 1516 . On Earth, storms frequently generate gust fronts above the surface. They occur 
when downdrafts, generally cooled by rain evaporation, reach the ground below the storm, producing 
cold surface currents, also called cold pools 18 . The gust front is the leading edge of the cold pool, and 
its propagation is mostly controlled by the direction of upper winds. 

To investigate whether Titan’s storms might lead to the production of such gusts with a prefer¬ 
ential direction, we analysed their formation under the wind-shear produced by Titan’s superrotation 
with a mesoscale regional model. We ran 2-dimensional (longitude-altitude) mesoscale simulations 
with TRAMS® 7 . This model includes full cloud microphysics (see Supplementary Information) with 
condensation, melting, freezing and evaporation of methane cloud particles. 

For the initial conditions, we used the temperature profile measured by Huygens — and a wind pro¬ 
file from the Titan IPSL GCM with no methane cycle 8 . Huygens’ measurements indicated a methane 
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Figure 1: 2D (altitude/longitude) simulation of the evolution of a storm under Titan’s conditions 
at equator during equinox, (a), (b), (c) and (d) are taken lh, lh40, 10h35 and 12h30 after the start of 
the simulation, respectively. The initial wind profile was derived from the GCM. The initial methane 
humidity corresponds to a CAPE of 500 J/kg. Color bars indicate the mixing ratio of condensed 
methane in g/kg. The wind vectors are scaled to the axis. A reference vector of 10 m/s for zonal wind 
is shown. In (c) and (d), the vertical wind component is reduced by 80 % to better see the gust front. 


humidity of 40% at the surface 20 . These conditions do not allow the development of convective 
clouds but were measured at the end of the southern summer. For our simulations at the equinoctial 
season, we used a methane profile with a surface humidity of ~50% or 60%, which is necessary to 
form convective clouds and corresponds to a convective available potential energy (CAPE) of 250 or 
500 J/kg. The deep convection was triggered by a warm bubble of +2 K (i.e. the air temperature is 
increased in the first km). In our simulations, convective clouds reach an altitude of 25 km where flow 
fast eastward winds, and the storm dynamics leads to the formation of a eastward gust front above the 
surface. Fig. |T] shows the development of a gust front in one of our simulations. The triggering of 
the convection is associated with symmetric convergence of air below the cloud (Fig. |T}.i). The cloud 
develops obliquely because of the wind-shear (Fig. El). To compensate for the ascending air in the 
storm, a downdraft forms behind the storm (Fig. ID>). Passing below the cloud, it is cooled by rain 
evaporation, which accelerates it, and flows at the surface, producing the gust front. Subsequently, 
the gust front ends up beyond the cloud (Fig. |T]j). It raises surrounding air from the surface to 3-4 
km, producing secondary convective clouds that can increase its energy and life-time (Fig. |Tj.i). The 
passage of the gust front also produces a local peak of eastward surface wind (Fig. |2ji), traveling over 
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large distances (in the order of 500 km in our simulations). Behind the gust front, wind speeds are 
typically one order of magnitude higher than usual winds close to the surface, and can reach 10 m/s 
during up to 9 hours (Fig. [2jr and b). 



Wind speed (m/s) 



Figure 2: Wind speed in the gust front, (a) Wind profiles in the gust front (red) and without storm 
(black dashed line), (b) Evolution of the zonal wind at 40 m (black line) during the passage of a 
storm. Eastward winds are positive. The red (blue) dashed line corresponds to the threshold friction 
speed of 0.04 m/s, i.e. 0.9 m/s at 40 m, for eastward (westward) 
winds. 

Our 2D simulations are representative of large storms or mesoscale convective systems. Individ¬ 
ual small storms (e.g. less than 30 km in diameter) should produce weaker and more isotropic gusts 
than in our simulations. Yet, multiple small storms should interact with each other and organize the 
convection to form large storm or mesoscale convective systems, as the giant cloud systems observed 
in September and October 2010 15 , producing strong gust fronts with an eastward propagation. The 
dynamics of such large storms should be essentially 2D and qualitatively well represented by our ide¬ 
alized 2D simulations (see Supplementary Information). We estimate that each region of the equator 
might statistically undergo on the order of one gust front every equinox (i.e. 2 per Titan year, see 
Supplementary Information). 

The nature of the sand on Titan is not well known. The optimum grain size for saltation is es¬ 
timated to be around 300 /tm with a threshold friction speed around 0.04 m/s 21 , corresponding to 
a wind speed around 0.9 m/s at 35 m above the surface (see Supplementary Information). In our 
GCM, winds at 35 m and at the equator are generally below 0.5 m/s and do not exceed 1 m/s (see 
Supplementary Fig. 6 and Fig. 7). Thus, usual winds on Titan are likely to generate only limited 
sand transport. Our GCM winds follow a Weibull distribution. They exceed the threshold friction 
speed of 0.04 m/s only 0.06% of the time, producing a westward sand flux at the equator of only 
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0.0015 m 2 /yr (see Supplementary Information for the sand flux calculation). However, GCMs tend 
to underestimate sand flux as they miss local gusts® 3 ! To compensate for this, we also consider a 
10 times stronger westward sand transport (0.018 m 2 /yr) by increasing GCM wind speed by 20%. 
Storm gusts largely exceed the wind threshold by an order of magnitude. According to our mesoscale 
simulations, the sand flux produced by storms is 2.7 times stronger eastward than westward, and one 
storm per equinox (i.e. two per Titan year) produces on average an eastward sand flux around 0.15 
m 2 /yr. Although episodic, tropical storms should therefore dominate the sand transport. A higher 
saltation threshold friction speed (e.g. 0.052 m/s) 22 would strongly reduce the sand flux from general 
circulation winds, while weakly impacting sand transport by storm gusts. 

On Earth, while in unlimited sediment many linear dunes are oriented orthogonally to the direction 
of maximum sediment transport® 4 , in limited sediment they align with resultant transport direction, 
called the resultant drift direction (RDD), and propagate over long distances with regular spacing 
and a limited amount of sediment in interdune areas- 2 ^. This mechanism may explain the formation 
and the evolution of Titan’s dunefields. It requires neither cohesion 26 nor complex secondary flows, 
but relies on the presence of sedimentary reservoirs and multidirectional wind regimes producing an 
overall sediment flux over a non-erodible ground 25 . Given that many locations in Titan’s dune fields 
appear to have sand-free bedrock interdunes 3 , the conditions of linear dune down-axis extension seem 
to be met. 

Fig. [3jt shows the RDD combining winds from the Titan IPSL GCM with no topography and 
storm gusts from mesoscale simulations, with respect to the storm frequency. With no storm or 
a very small frequency of storm, the RDD is oriented westward with a southward component due 
to Saturn’s eccentricity (the latter results in stronger solar forcing and increased Hadley circulation 
during the southern summer, hence producing southward winds stronger than northward winds during 
the northern summer). With around two storm per Titan year (i.e. one per equinox), the sand tranport 
produced by equinoctial storms dominates and the RDD is oriented eastward. With such a frequency, 
the storm control remains efficient even for the GCM winds with increased gust. 

For a 100 m-high dune, the eastward sand flux produced by two storms per Titan year (i.e. 0.15 
m 2 /yr) corresponds to an extension rate of 3 mm/yr (see Supplementary Information). According to 
the typical length of Titan’s dunes (i.e. 30-50 km) 27 , this implies a minimum formation time of around 
15 Myr. This time is longer than the period of Saturn’s perihelion precession (i.e. 45 kyr), which 
corresponds to the switch of the warmest summer between both hemispheres. Thus the southward 
component due to this orbital effect should vanish in dune orientation. Barchans and small dunes 
have been observed with shifted orientation as compared to long dunes 28 . This has been interpreted 
as the effect of a long-term, cyclic multimodal wind regime. While small dunes may respond to long 
orbital changes (affecting both general circulation winds and storm impact), our multimodal wind 
pattern can also explain their shifted orientation by growth orthogonal to the direction of maximum 
sediment transport 24 25 . Fig. |3j? and Fig. [3]e show the RDD map in the equatorial band without or with 
the effect of two storms per Titan year. In both cases, we averaged the effect of Saturn’s eccentricity 
according to the above discussion. Our GCM predicts a RDD slightly converging to the equator 
for latitudes lower than 12° N/S and diverging from the equator for latitudes higher than 12° N/S. 
The divergence is due to faster poleward summer winds. This may explain the observed poleward 
orientation of dunes for latitudes higher than around 10° N/S 229 . These poleward sand fluxes also 
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Figure 3: Storm impact on the resultant drift direction (RDD). (a) The RDD as a function of the 
storm frequency per equinox with a saltation threshold of 0.04 m/s, using GCM winds (black solid 
line) or GCM winds with increased gust (black dashed line). Angles are measured anti-clockwise 
from the East (i.e. -180° is westward, -90° is southward and 0° is eastward). The red line symbol¬ 
izes the passage from westward to eastward dune growth, (b) Map dune orientation observed with 
Cassini’s radar® 9 (red) and map of RDD obtained with the GCM (with increased gust) for low lati¬ 
tudes with a threshold of 0.04 m/s and no storm effect (black), (c) is same as (b) but with the impact 
of two storms per Titan year. 


imply that Titan’s sand cannot be transported from polar regions to the equatorial band. The sand of 
Titan’s dunes must therefore have been produced in equatorial regions. 

The sand flux pattern that we estimate for Titan’s dunes is similar to those implicated in the 
formation of longitudinal dunes in the Rub’al-Kali desert (southern Arabia) and in the Great Sand 
Sea (Egypt). In these deserts subject to a dominant wind, linear dunes are parallel to the RDD (see 
Fig. [4j> and Supplementary Fig. 10). We can therefore establish an analogy with Titan’s dune 
formation (see Fig. [4]). In terrestrial deserts, gust fronts formed by convective clouds can produce 
dust storms, also called haboobs 30 . Similarly, gust fronts generated by Titan’s storms could produce 
dust storms. Indeed, their friction speed might exceed the threshold for micrometric particles (around 
4 m/s at 35 m for a 10 /im diameter particles, see Supplementary Information). They could also 
transport centimetric pebbles (threshold speed around 4-9 m/s at 35 m for 1-5 cm diameter pebbles), 
explaining the diffuse material trailing out to the east around many mountains 27 . 

In conclusion, the combination of GCM and mesoscale simulations, together with analogies with 
terrestrial dune formation, provides a comprehensive scheme to explain the major features of Titan’s 
dunes (i.e. eastward propagation, linear shape, poleward divergence and 2-3 km spacing). This anal¬ 
ysis highlights the potential importance of storms and mesoscale phenomena on aeolian processes. 
Such rare and strong events may play or have played a important role in some dune fields on Earth 
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Figure 4: Analogy between linear dunes on Titan and in Rub’al-Kali desert on Earth, (a) de- 

noised image of Titan’s dunes from Cassini’s radar SAR (see Supplementary Information). The inset 
shows the sand flux rose similar to Fig. |3jt, but calculated at 7.5°S and averaging the effect of Sat¬ 
urn’s eccentricity, (b) longitudinal dunes in Rub’al-Kali desert (18° N, 48° E) with the sand flux roses 
calculated from winds at 10 m. For both images, the arrow corresponds to the resultant drift direction 
calculated with the sand flux rose. 

and Mars. 

Methods 

Our analysis was performed using GCM and mesoscale simulations, together with analogies with 
terrestrial dunes. We use the Titan IPSL GCM- 10 to predict surface winds in the equatorial regions. 
This GCM reproduces well the superrotation and the thermal structure measured by the Huygens 
probe in the troposphere. Simulations are performed taking into account the effects of the diurnal 
cycle and Saturn’s gravitational tides, but with no topography and no methane cycle (see also Sup¬ 
plementary Information). Titan’s climate is essentially dry, excepted in polar regions. The methane 
cycle is therefore not expected to significantly affect the tropospheric dynamics in the equatorial re¬ 
gions. This is illustrated by the excellent match between the thermal structure and its different layers 
measured by the Huygens probe and simulated by the GCM®. To compute wind statistics and sand 
transport, we use instantaneous GCM surface winds over a year and all longitudes. 

Methane storms are simulated with the mesoscale model TRAMS in 2D including full methane 
cloud microphysics (see also Supplementary Information) 67 . Incidentally, a GCM including the 
methane cycle cannot resolve storm cold pools and gust fronts, which are subgrid phenomena. As 
boundary conditions, we use the GCM zonal wind profile at the equator and at the north spring 
equinox, the Huygens’ temperature profile 19 and the Huygens’ methane humidity profile with an in¬ 
crease in the lower troposphere. These idealized simulations are representative of large or squall line 
storms. We estimate the order of magnitude of the storm frequency at each location in the equatorial 
band around one storm every equinox. A simple estimation of the size and displacement of the giant 
arrow shape storm observed by Cassini at the equinox in 2009^ 1516 gives a lower limit of 0.2 storm 
per equinox (see Supplementary Information). To compute the sand transport and sand flux roses 
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produced by general circulation winds and gust fronts, we use a sand flux formula with a saltation 
threshold speed estimated from previous studies 21 . We also establish an analogy between Titan’s 
dunes and some terrestrial dune fields, in particular in southern Arabia and Egypt (see Fig. [4] and 
Supplementary Fig. 10), where we predict sand fluxes from wind reanalysis data. The simulation 
data are available from B.C. on request. 
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Methane storms as a driver of Titan’s dune orientation 

Supplementary Information 

1 Orientation of Titan’s dunes from radar images 



Figure 5: Map of Titan’s dune orientation. Map of radar-measured dune orientation vectors 2 ' 
(CREDITS: NAS A/J PL-Cal tech/A S I/S pace Science Institute, PIA11801), showing the global east¬ 
ward propagation and the divergence from the equator for latitudes higher than 10°. 

Fig. [5] shows the map of dune orientation and propagation, indicated with vectors, over a near- 
infrared basemap derived from Cassini ISS (Imaging Science Subsystem) images 2 . The direction 
of propagation is obtained by looking at dune morphology around obstacles and dune termina¬ 
tions, with for instance some dunes stopping on the west side of obstacles and others diverting 
and recombining beyond the east side of obstacles. This analysis reveales that all dunes propagate 
eastward. Concerning the North-South orientation, a statistical analysis indicates that dunes tend 
to diverge poleward for latitudes higher than around 10° N/S 29 . 

Images of Titan’s surface obtained by the Cassini’s radar SAR (Synthetic Aperture Radar) suffer 
from errors from a variety of sources, the most prominent being speckle noise. This noise refers 
to the constructive and destructive interference of scattered energy from roughness elements on 
a scale smaller than the size of a SAR pixel. The result is a multiplicative noise, which hinders 
interpretation. An advanced denoising algorithm has been adapted to Cassini SAR data 31 . It has 
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been applied to the original image (i.e. the image by Cassini Radar, T8 flyby) of Fig. 4a of the 
letter to significantly reduce the noise. 

Fig. 3a and 3b of the article show dune orientations obtained in the same way as Fig. [5]but with 
denoised radar images whose dunes segments have been detected 29 and their orientation averaged 
over areas of 8° x 8° in longitude-latitude. 


2 Calculation of general circulation winds 


Description of the IPSL Titan GCM 

For this study, we used a 3-dimensional GCM 810 . A horizontal resolution of 32 x 48, correspond¬ 
ing to resolutions of 3.75° latitude by 11.25° longitude, is used for the simulations. This GCM 
covers altitudes from the ground (first level at 35 m) to 500 km. The dynamical core is based on 
the most up-to-date version of the LMDZ 32 . It is a finite-difference discretization scheme that 
conserves both potential enstrophy for barotropic nondivergent flows, and total angular momen¬ 
tum for axisymmetric flows. The version used in this study includes gravitational tides 33 , though 
the impact in the troposphere does not influence the effects described in this work. We found tidal 
effects on the pressure similar to previous works^H but tidal winds are much weaker in our model. 
We use a fully coupled aerosol microphysics calculated in 2D (zonally averaged) 34 . The present 
model is dry and does not take into account the methane cycle. The profile of methane is fixed 
(close to the HASI profile 2 ®) for the radiative transfer. The latter is based on the McKay radiative 
code 35 and includes the diurnal cycle. The radiative transfer is called 200 times per Titan day. For 
the surface, we use an albedo of 0.3, a rugosity length of 0.005 m, an emissivity of 0.95 and a 
thermal inertia of 400 J/m 2 /K. In this study, we run the GCM with a flat topography. 

To calculate the friction speed u* from the GCM wind u at an altitude z, we use the relation: 


it* 


K 


In ( z / zq ) 


u 


( 1 ) 


with k = 0.4 the Von-Karman constant and zq = 0.005 m the rugosity length. The threshold 
u* t = 0.04 m/s corresponds to a wind speed of 0.89 m/s at 35 m. 


GCM wind statistics: 

For all calculation implying the GCM wind statistics (e.g. sand fluxes), we use instantaneous 
GCM winds with 20 outputs per Titan day and combining wind series at all longitudes. The 
dissimilarities for wind statistics obtained for different Titan years are negligeable. Fig. [6] 
shows the surface wind roses produced by the GCM in the equatorial band. Surface winds are 
essentially bimodal. At the equator, they blow from NE to SW in northern winter and from SE 
to NW in northern summer. For latitudes higher than 5° N/S, the summer winds are eastward. 
Because of Saturn’s eccentricity, the southern summer is shorter and hotter than the northern 
summer. This implies that northerly winds are less frequent but stronger than southerly winds. 
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Figure 6: Wind roses (direction and speed) produced by the GCM. The roses have been obtained 
with instantaneous winds at 35 m at 0°N (left) 10°N (middle) and 20°N (right). The frequency is given 
in percent. 
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Figure 7: Statistics of GCM winds at the equator. The different lines correspond to the relative 
probability of friction wind speed per bin of 0.001 m/s for the GCM winds (in red) and the GCM winds 
with increased gust (in orange) at the equator. The black line corresponds to the Weibull distribution 
fitting the GCM wind statistics (coefficient k=2.15 and c=0.0157 m/s) and the blue line the distribution 
fitting the GCM wind statistics with increased gust (coefficient k=2.15 and c=0.019 m/s). 


Fig. [7] shows the relative probability of the friction speed from the GCM (in red) at the equator 
and per bin of 0.001 m/s. Wind speed exceeds the threshold friction speed of 0.04 m/s only 
around 0.06 % of the time. The GCM wind statistics are well described by a Weibull distribution 
(in black in the figure), for which the probability of exceeding a friction speed U is P(> U) = 
exp(— (U/c) k ), with c=0.0157 m/s the scale parameter and k= 2.15 the shape parameter. In order 
to represent the higher variability of Titan’s winds due to local gusts that are not captured with 
the low-resolution GCM grid (typical size of the spatial grid: 500 kmx 170 km), we have also 
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considered the case of an increase of wind speed by 20% (in orange in the figure). This arbitrary 
value corresponds to a significant increase of wind and leads to a one order of magnitude stronger 
sand transport. Because of the weak turbulence in the boundary layer®, it is likely to be a quite 
high value for high value for the representation of missing gusts. Concerning wind statistics, it 
is equivalent to an increase of the scale parameter of the Weibull distribution by 20% (in blue in 
the figure). In these conditions, the wind speed exceeds the threshold of 0.04 m/s around 0.7 % 
of the time (i.e. around ten times more than without correction). 


3 Simulation of methane storms 

Description of the TRAMS model 

The Titan Regional Atmospheric Modeling System (TRAMS) is a coupled, regional-scale dy¬ 
namics and column microphysics model® 7 . The governing equations for the dynamical core are 
the standard non-hydro static Reynolds-averaged primitive equations 37 . The microphysics pack¬ 
age was adapted from Barth et al. (2006) 38 and operates independently on each column. Methane 
cloud particles form through nucleation onto submicron-sized haze particles. Through conden¬ 
sation and coalescence, methane particles can grow up to millimeter sizes. Depending on the 
temperature of the surrounding environment, the methane clouds form as either ice particle or 
droplets; melting and freezing of cloud particles is also included. Liquid droplets are treated as a 
mixture of N2/CH4 following the work of Thompson et al. (1992) 39 . 

TRAMS is run here as a 2-D model. The 2-D simulations are run with a horizontal grid spacing 
of 1 km, and a total horizontal extent of 1000 km. Cyclic horizontal boundary conditions are 
employed. The vertical grid extends up to about 50 km; vertical layers are more compact near the 
surface, starting at about 15 m spacing and extend to constant 2 km spacing above 15 km altitude. 
The atmosphere is initialized horizontally homogeneously using the temperature-pressure profile 
measured by the HASI instrument on the Huygens probe 19 . An initial horizontal wind is 
included using the w-velocity component from the GCM simulations. For methane, we construct 
profiles using a fixed amount of convective available potential energy (CAPE); we look at 
cases with CAPE=250 and CAPE=500 (equivalent to a methane mixing ratio of 5 g/kg, or 10 
g/kg, respectively, near the surface). Cloud formation is initially triggered by perturbing the 
atmosphere with a warm bubble (the air is locally warmed in the first km), which has a maximum 
atmospheric temperature increase of 2 K. This 2K perturbation is a classical method to trigger 
deep convection in mesoscale simulations of terrestrial storm! 40 ! It is also possible to trigger 
deep convection with initial vertical updrafts (typically 1 m/s) but both methods lead to similar 
cloud dynamics 41 . In reality, the moist convection would be triggered by updrafts produced by 
planetary waves, the Hadley cell convergence, the diurnal cycle or the topography. The updraft 
velocities at largescale trigerring storms could be obtained from a GCM. 

Dynamics of convective cloud systems 

2D simulations of methane storms have been performed with our mesoscale model investigating 
in more details the effect of wind shear and CAPE on the storm dynamics, morphology and 
life time 42 . It revailed behaviors similar to 2D simulations of terrestrial storms 43 . We therefore 
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expect the dynamics of the various kinds of Titan’s storms (i.e. from individual convective cell 
to mesoscale convective systems) to be very similar to those of terrestrial storms. In particular, 
in our simulations the propagation and lifetime of Titan’s storms is mostly controlled by the gust 
front. The leading edge of the gust front is able to lift moist air from the surface to trigger a 
new cell. When this new cell is close enough to the previous cell they merge together, increasing 
the lifetime of the storm and the strength of the gust front. When the gust front moves too 
quickly compared to the storm, they become disconnected and the storm dissipates. If a new 
cell were produced it would not be able to merge with the previous one. Fig. 1 of the letter in 
the manuscript reveals this behaviour. A new cell is produced in Fig. lc at around 20 km of 
the pre-existing cell and merges with it. Other cells were produced by the leading edge of the 
gust front before this event. We found that these new cells form at 20-40 km in front of the 
pre-existing cell before merging with it. In figure 4d, the gust front has moved away from the 
mean cell. It triggers a new cell but too far (100 km) from the pre-existing cell and they dissipate. 
Our idealized 2D simulations of methane storms and gust fronts are representative of large 
storms or mesoscale convective systems (typically 100-1000 km in latitude), as those observed 
by Cassini 1516 . Small isolated convective clouds, as the ones simulated in 3D by Hueso and 
Sanchez-Lavega 2006 41 , should produce weaker gust fronts than in our simulations and with 
a more istropic direction. However, multiple small convective clouds should merge into larger 
storms (see supplementary figure 2 in Hueso and Sanchez-Lavega 2006) associated with stronger 
gust fronts and with a spanwise dimension similar to the one of these cloud systems (i.e. 100-1000 
km). We also expect that some large storms or mesoscale convective systems on Titan evolve to 
produce squall lines or bow echos as they generally do on Earth®®. 

On Earth, a mesoscale convective system has a 3D structure, yet it becomes mostly 2D when 
the wind shear is essentially 2D 18 40 45 46 . On Titan, the meridional wind shear is very weak com¬ 
pared to the longitudinal wind shear corresponding to the superrotation. Moreover, the Coriolis 
force is particularly weak at low latitude on Titan. This implies that equatorial mesocale con¬ 
vective systems should propagate westerly, forming fronts with no deformation by the Coriolis 
force. We therefore expect Titan’s storm dynamics to be essentially 2D. This implies that surface 
winds produced by Titan’s storms are primarily oriented east-west and that gust fronts should be 
qualitatively well represented by our idealized 2D simulations. 

Finally, the direction of propagation of the gust front produced by a mesoscale convective system 
is essentially controlled by the transport of momentum from the mid-troposphere to the surface 
by cold pools^ 8 . The presence of the superrotation necessarily implies a favoured eastward propa¬ 
gation for gust fronts on Titan. A good analogy with the generation of eastward surface winds by 
Titan’s storms is the impact of mesoscale convective systems over the Pacific warm pool during 
the Madden-Julian Oscillation. During such periods, mesoscale convective systems accerelate 
surface winds eastward in regions of fast high eastward winds 47 49 . 


4 Estimation of storm frequency 

If we consider only the giant storm which was observed on 27 September 2009 1 ' and if we 
assume that the same event occurs every equinox, we can evaluate a lower limit for the storm 
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frequency as: 


9 A x 

^^T-storm 

J min ~ ~1 \A) 

-^•equator -L Titan 

A storm is the area covered by the passage of the storm, A equat or ~ 4x 10 7 km 2 is the area of the 
equatorial band (30° S to 30° N) and T Titan is the length of a Titan year (~29.5 terrestrial years). 
The width of the storm was around 2000 km and it traveled at least 4000 km during probably a few 
days to produce the observed surface changes 16 . Thus, thee passage of this storm covered around 
20% of Titan’s equatorial band. This gives an average storm frequency of 0.4 per Titan’s year or 
0.2 storm per equinox. It seems reasonable to consider this value as a lower limit for the storm 
frequency. Other large cloud systems have been observed during this equinox 12 ®^ 17 50 . Moreover, 
Titan has only been observed during a short period of the equinox and not globally. Possible 
observations of liquid ethane on Titan’s surface suggest that moderate and small rainstorms are 
quite common during the equinoctial season at the equator 51 . We therefore expect that the mean 
storm frequency should rather be in the order of one storm per equinox (i.e. 2 storms per Titan 
year or 0.068 storms per terrestrial year). 


5 Sand flux calculations 


Saltation threshold The saltation threshold corresponds to the minimal friction speed for which 
the wind stress is sufficient to lift particle 54 . It has been estimated to be around 0.04 m/s for 
Titan 21 52 54 for an optimum particle diameter of 200-300 pm and a sediment density of around 
1000 kg/m 2 . A simple expression of the threshold friction speed for saltation is^ 
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with Atv=0.111 a dimensionless parameter, p air the air density, p se< j the sediment density, I) the 
particle diameter, g Titan’s gravity, 7 a parameter scaling the strength of the interparticle forces. 
An actually depends on the particle friction Reynolds number 5 ^_but this dependence remains 
limited and An can be considered constant at first approximation 55 . 

Fig. [ 8 ] shows the threshold friction speed according to particle diameter for g = 1.35 m s -2 , 
p se d=1000 kg/m 3 , p air =5.3 kg/m 3 and 7 = 1 .5 x 10 _ 4 N/m. For these values, the minimal threshold 
friction speed is 0.045 m/s for an optimum diameter of 350 pm. For pebbles with a diameter of 
1-5 cm, the threshold friction speed is 0.18-0.4 m/s. This corresponds to a wind speed at 35 m of 
4-9 m/s (see below). For dust particles with a diameter of 10 /mi, the threshold friction speed is 
0.19 m/s. This corresponds to a wind speed at 35 m of around 4 m/s. General circulation winds 
never reach such high thresholds. In contrast, gust fronts produced by methane storms may be 
strong enough to transport centimetric pebbles and dust particles. Dust particles can also be raised 
by the impact of saltating sand, thus with a wind speed lower than threshold friction speed. The 
release of dust particles could lead to the formation of dust storms that could persist in the lower 
atmosphere several hours or days after the passage of the methane storm. 


16 











Figure 8 : Threshold friction speed for saltation as a function of the particle diameter. The salta¬ 
tion speed has been calculated using relation (3) for Titan’s conditions 2 ', with a sediment density of 
1000 kg/m 3 and a parameter y-1. 6 x 10 _ 4 N/m. 

Sand transport law 

For calculating sand transport, we use the law from Kawamura (1951) 53 57 58 : 

Q = 2.6 (J (w* - + u * t ) 2 (4) 

\ 9 Psed J 

with Q the sand flux per unit of width in m 2 /s, p :ur = 5.3 kg/m 3 the air density, p sed = 1000 kg/m 3 
the sediment density, g = 1.35 m s ' 2 Titan’s gravity, w* the friction speed and u xi the threshold 
friction speed for transport. 

To calculate //„ from u, the GCM or meso-scale wind at an altitude z, we use the relation (1). 
Several other sand transport laws have been proposed^ 54 ! Their uses instead of the Kawamura’s 
law could quantitively change our results (e.g. mean sand flux values) but not qualitatively (e.g. 
the reorientation of sand flux by methane storms). 

Calculation of the storm impact on the sand flux 

To calculate the impact of storms on the sand flux, we integrate formula (4) for the sand flux 
produced by the passage of one storm in our simulation and averaged it over all the domain 
(1000 km). In our simulations, storms are around 40 km width (in longitude) and travel around 
200 km during around one day. Thus, we multiply the previous value by 5 (1000 km/200 km) to 
get the average impact of one storm. Finally, we divide this value by a period of half a Titan year. 
This yields to the mean sand flux produced by one storm per equinox (typically 0.15 m 3 /m/year). 
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Figure 9: Sand flux roses obtained by combining the GCM winds with winds produced during 
one typical gust front every equinox. The arrows correspond to the resultant drift directions, a, b 
and c correspond to 0°, 10° and 20° N latitudes, respectively, with GCM winds (speed increased by 
20 %) and a threshold of 0.04 m/s, and averaging the effect of Saturn’s eccentricity. 


To obtain the total sand flux rose, we added the eastward and westward component of this sand 
flux, multiplied by the storm frequency per equinox, to the rose obtained with the GCM. In our 
mesoscale simulations, the eastward sand flux is around 2.7 times higher than the westward sand 
flux. The direction of surface winds produced by real storms is not likely to be purely west-east 
as in our 2D simulations. However, the favoured direction is eastward. Thus, the sand flux and 
the RDD are eastward on average. To get more realistic figures, we added a small ad hoc angular 
dispersion, for the sand flux rose produced by storms: Q oc e and 0 = 0 for 0 |>20° where 

6 is the direction. Fig. [9] shows to the sand flux roses for different latitudes, averaging the effect 
of Saturn’s eccentricity. 

Longitudinal dune extension rate The extension rate of a longitudinal dune can be expressed 
as: 

W 

c = Qrbp x (5) 

c the extension rate (in m/s), Qrdp is the norm of the mean flux, also called the resultant drift 
potential (in m 2 /s), W is the width of the dune (in m) and S is the cross section (in m 2 ). For a 
linear dune, S — W x H /2 (i.e. the surface of an isocele triangle) where H is the height of the 
dune. We have then c = 2Q KDP /H. 

For Titan’s dunes, H m 100 m. If we consider the impact of one storm every equinox, Qrdp = 
0.15 m 2 /yr and the extension rate is around 3 mm/yr. 

Sand flux in southern Arabia and in Egypt To calculate the sand flux in Rub’al-Kali 
desert (Fig. 4b of the letter) and in the Great Sand Sea (Fig. [TO]), we use the wind outputs of 
the ERA-Interim-project, a reanalysis from the European Centre for Medium-Range Weather 
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Forecasts 159 * We analyzed the 10 m wind data for the period between the 1/1/1979 and the 
31/12/2012, at 18°N and 48°E for Rub’al-Kali desert, and at 25.5°N and 26.25°E for the Great 
Sand Sea. Mean sand flux roses (see Fig. 4b of the letter and Fig. [TO]) have been produced using 
these wind data and the formula (4). The mean sand flux direction (i.e. the RDD) is shown with 
an arrow and is parallel to the longitudinal dune orientation as predicted* 25 } 



Figure 10: Longitudinal dunes in Egypt and Algeria. Map data: Google, (a) Dune field in Egypt 
(25.5°N, 26.25°E) with the sand flux roses calculated from winds at 10 m. (b) Dune field in Algeria 
(29.5° N, 5.5° E). The arrow corresponds to the resultant drift direction calculated with the sand flux 
rose. 
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